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Planetary turbulence is observed to self-organize into large-scale structures such 
as zonal jets and coherent vortices. One of the simplest models that retains the 
relevant dynamics of turbulent self-organization is a barotropic flow in a beta-plane 
channel with turbulence sustained by random stirring. Non-linear integrations 
of this model show that as the energy input rate of the forcing is increased, the 
homogeneity of the flow is first broken by the emergence of non-zonal, coherent, 
westward propagating structures and at larger energy input rates by the emergence 
of zonal jets. The emergence of both non-zonal coherent structures and zonal jets 
is studied using a statistical theory. Stochastic Structural Stability Theory (S3T). 
S3T directly models a second-order approximation to the statistical mean turbulent 
state and allows the identification of statistical turbulent equilibria and study of 
their stability. Using S3T, the bifurcation properties of the homogeneous state in 
barotropic beta-plane turbulence are determined. Analytic expressions for the zonal 
and non-zonal large-scale coherent flows that emerge as a result of structural instability 
are obtained and the equilibration of the incipient instabilities is studied through 
numerical integrations of the S3T dynamical system. The dynamics underlying 
the formation of zonal jets are also investigated. It is shown that zonal jets form 
from the upgradient momentum fluxes that result from the shearing of the eddies 
by the emerging infinitesimal large-scale flow. Finally, numerical simulations of the 
nonlinear equations confirm the characteristics (scale, amplitude and phase speed) of 
the structures predicted by S3T, even in highly non-linear parameter regimes such as 
the regime of zonostrophic turbulence. 


Atmospheric and oceanic turbulence is commonly ob¬ 
served to be organized into spatially and temporally co¬ 
herent structures such as zonal jets and coherent vortices. 
A simple model that retains the relevant dynamics, is a 
barotropic flow on a /3-plane with turbulence sustained by 
random stirring. Numerical simulations of the stochasti¬ 
cally forced barotropic vorticity equation on the surface of a 
rotating sphere or on a /3-plane, have shown the coexistence 
of robust zonal jets and of large-scale westward propagating 
coherent structures that are referred to as satellite modes 
(Danilov and Gurarie 2004) or zonons (Sukariansky et al. 
2008; Galperin et al. 2010). Emergence of these coherent 
structures in barotropic turbulence has also another feature. 
As the energy input of the stochastic forcing is increased or 
dissipation is decreased, there is a sudden onset of coher¬ 
ent zonal flows (Srinivasan and Young 2012; Constantinou 


et al. 2014) and non-zonal coherent structures (Bakas and 
loannou 2014). This argues that the emergence of coherent 
structures in a homogeneous background of turbulence is a 
bifurcation phenomenon. 

An advantageous method to study such a phenomenon, 
is to adopt the perspective of statistical state dynamics of 
the flow, rather than look into the dynamics of sample real¬ 
izations of direct numerical simulations. This amounts to 
study the dynamics and stability of the statistical equilibria 
arising in the turbulent flow, which are fixed points of the 
equations governing the evolution of the flow statistics. This 
approach is followed in the Stochastic Structural Stability 
Theory (S3T) (Farrell and loannou 2003) or Second Order 
Cumulant Expansion theory (CE2) (Marston et al. 2008). 
This theory is based on two building blocks. The first is 
to do a Reynolds decomposition of the dynamical variables 
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into the sum of a mean value that represents the coherent 
flow and fluctuations that represent the turbulent eddies 
and then form the cumulants containing the information on 
the mean values (first cumulant) and on the eddy statistics 
(higher order cumulants). The second building block is 
to truncate the equations governing the evolution of the 
cumulants at second order by either parameterizing the 
terms involving the third cumulant (Farrell and loannou 
1993a,b,c; DelSole and Farrell 1996; DelSole 2004) or setting 
the third cumulant to zero (Marston et al. 2008; Tobias 
et al. 2011; Srinivasan and Young 2012). Restriction of 
the dynamics to the first two cumulants is equivalent to 
neglecting the eddy-eddy interactions in the fully non-linear 
dynamics and retaining only the interaction between the 
eddies with the instantaneous mean flow. While such a 
second order closure might seem crude at first sight, there 
is strong evidence to support it (Bouchet et al. 2013). 

Previous studies employing S3T have already addressed 
the bifurcation from a homogeneous turbulent regime to 
a jet forming regime in barotropic /3-plane turbulence and 
identified the emerging jet structures both numerically (Far¬ 
rell and loannou 2007) and analytically (Bakas and loannou 
2011; Srinivasan and Young 2012) as linearly unstable modes 
to the homogeneous turbulent state equilibrium. It was also 
shown that the resulting dynamical system for the evolution 
of the first two cumulants linearized around the homoge¬ 
neous equilibrium possesses the mathematical structure of 
the dynamical system of pattern formation (Parker and 
Krommes 2013). Comparison of the results of the stability 
analysis with direct numerical simulations have shown that 
the structure of zonal flows that emerge in the non-linear 
simulations can be predicted by S3T (Srinivasan and Young 
2012; Constantinou et al. 2014). However, these research 
efforts, have assumed that the ensemble average employed 
in S3T is equivalent to a zonal average, a simplification that 
treats the non-zonal structures as incoherent and cannot 
address their emergence and effect on the jet dynamics. In 
addition, the eddy-mean flow dynamics underlying the S3T 
instability even in the jet formation case, that involve only 
the interactions of small scale waves with the large-scale 
coherent structures are not clear. 

So the goals in this article are the following. The first 
goal is to adopt a more general interpretation of the ensem¬ 
ble average, in order to address the emergence of coherent 
non-zonal structures. We adopt the more general inter¬ 
pretation that the ensemble average is a Reynolds average 
over the fast turbulent motions (Bernstein 2009; Bernstein 
and Farrell 2010). With this definition of the ensemble 
mean, we obtain the statistical dynamics of the interaction 
of the coarse-grained ensemble average field, which can 
be zonal or non-zonal coherent structures represented by 
their vorticity, with the fine-grained incoherent field repre¬ 
sented by the vorticity second cumulant and we revisit the 
structural stability of the homogeneous equilibrium under 


this assumption. The second goal is to study in detail the 
eddy-mean flow dynamics underlying the S3T instability 
focusing on the example of jet formation. And the third 
goal is to compare the characteristics of the structures that 
emerge in S3T against non-linear simulations, even in highly 
non-linear regimes that at first glance present a challenging 
test for the restricted dynamics of S3T. 

1. Formulation of Stochastic Structural Stability 

Theory under a generalized average 

Consider a nondivergent barotropic flow on a /3-plane 
with cartesian coordinates x = {x,y). The velocity field, 
u = (u, v), is given by (m, v) = {—dytp, 9x4’)^ where ip is the 
streamfunction. Relative vorticity (P{x,y,t) = Aip, evolves 
according to the non-linear (NL) equation: 

{dt + u-V)C + Pv = -rC-iyA'^C + y/ef‘', ( 1 ) 

where A = + dyy is the horizontal Laplacian, /3 is the 

gradient of planetary vorticity, r is the coefficient of linear 
dissipation that typically parameterizes Ekman drag in plan¬ 
etary atmospheres and v is the coefficient of hyper-diffusion 
that dissipates enstrophy flowing into unresolved scales. 
The exogenous forcing term /®, parameterizes processes 
such as small scale convection or baroclinic instability, that 
are missing from the barotropic dynamics and is necessary 
to sustain turbulence. We assume that is a temporally 
delta correlated and spatially homogeneous random stir¬ 
ring injecting energy at a rate e and having a two-point, 
two-time correlation function of the form: 

{r{xi,yi,ti)f{x2,y242)) = S{t2 -ti)E{xi,X2,yi,y2), 

( 2 ) 

where the brackets denote an ensemble average over the 
different realizations of the forcing. 

S3T describes the statistical dynamics of the first two 
same time cumulants of (1). The equations governing the 
evolution of the first two cumulants are obtained as follows. 
We decompose the vorticity field into the averaged field, 
Z = T[C], defined as a time average over an intermediate 
time scale and deviations from the mean or eddies, (p' = 
C — Z. The intermediate time scale is larger than the time 
scale of the turbulent motions but smaller than the time 
scale of the large scale motions. With this decomposition 
we rewrite (1) as: 

{dt + U ■W)Z + /3V = -W ■ r[u'C'] -rZ- (3) 

where U = \U,V] = [—dy'^^dx'^] and u' = [u',v'] = 
[—dyip', dxip'] are the mean and the eddy velocity fields 
respectively. The mean vorticity is therefore forced by the 
divergence of the mean vorticity fluxes. The eddy vorticity 
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(' evolves according to: 


where 


(a* + U • V) c' + (/3 + dyZ)v' + u'd^Z = 

= -rC - uA\' + f+ T[u' ■ VC'] - u' • VC', (4) 

'- ,, -^ 

where /"* is the term involving the non-linear interactions 
among the turbulent eddies. A closed system for the sta¬ 
tistical state dynamics is obtained by first neglecting the 
eddy-eddy term to obtain the quasi-linear system, 

{dt+V ■V)Z + (3V = -W ■ r[u'C'] - rZ - i^A^Z, (5) 

(a* + U • V) C + W + dyZ)v' + u'd^Z = 

= -rC - vA'^C' + Ve/®, (6) 


In order to obtain the statistical dynamics of the quasi- 
linear system (5)-(6) we adopt the general interpretation 
that the ensemble average over the forcing realizations 
is equal to the time average over the intermediate time 
scale (Bernstein 2009; Bernstein and Farrell 2010). Under 
this assumption, the slowly varying mean flow Z is also 
the first cumulant of the vorticity Z = (C), where the 
brackets denote the ensemble average. The time mean of 
the vorticity flux is equal to the ensemble mean of the 
flux T[u'C'] = (u'C'). The fluxes can be related to the 
second cumulant C(xi,X 2 ,t) = (C'(xi)C'(x 2 )), which is the 
correlation function of the eddy vorticity between the two 
points Xi = {xi,yi), i = 1,2. We hereafter indicate the 
dynamic variables that are functions of points x^ = {xi,yi) 
with the subscript i, that is C) = C'(^i)- By making the 
identification that the fluxes at point x are equal to the 
value of the two variable function (u^C^) evaluated at the 
same point x = xi = X 2 , we write the fluxes as: 


A, = -\J,-\/,-{l3 + dy^Z)d,^A-^ + d,^Zdy^A-^-r-i^Al 

( 11 ) 

governs the dynamics of linear perturbations about the in¬ 
stantaneous mean flow U. The right hand side of (10) is the 
correlation of the external forcing with vorticity, which for 
delta correlated stochastic forcing is independent of the state 
of the flow and is equal at all times to the prescribed forcing 
covariance: ^/e (ffC^ + /ICi) = e (/i/I) = Therefore 
The second cumulant evolves then according to: 

dtC = {Ai + A2)C + ( 12 ) 

and forms with Eq. (9) the closed autonomous system of 
S3T theory that determines the statistical dynamics of the 
flow approximated at second order. 

The S3T system has bounded solutions (cf. Appendix A) 
and the fixed points Z^ and , if they exist, define sta¬ 
tistical equilibria of the coherent structures with vorticity, 
Z^, in the presence of an eddy field with second order 
cumulant or covariance, . The structural stability of 
these statistical equilibria addresses the parameters in the 
physical system which can lead to abrupt reorganization of 
the turbulent flow. Specifically, when an equilibrium of the 
S3T equations becomes unstable as a physical parameter 
changes, the turbulent flow bifurcates to a different attrac¬ 
tor. In this work, we show that coherent structures emerge 
as unstable modes of the S3T system and equilibrate at 
finite amplitude. The predictions of S3T regarding the 
emergence and characteristics of the coherent structures 
are then compared to the non-linear simulations of the 
stochastically forced barotropic flow. 


(u'C') = (u'lC^).,.., . ( 7 ) 

Expressing the velocities in terms of the vorticity [u',v'] = 
\—dyA~^,dxA~^]C,', where is the integral operator that 
inverts vorticity into the streamfunction held, we obtain 
the vorticity huxes as a function of the second cumulant, 
in the following manner: 

(u'C') = [(WlC2)xi=X2 > ('^iC2)xi=X2] 

= - {a,. ,{a„ 

Consequently, the hrst cumulant evolves according to: 

dtZ + UZ,^ + U(/3 -f Zy) +rZ + vA^Z = 

= d. K - dy . (9) 

Multiplying (A5) for dtQ'i by C 2 and (A5) for dtC,^ t>y Ci, 
adding the two equations and taking the ensemble average 
yields the equation for the second cumulant (7: 

dtC - (Ai+ a2)c = (/rc^ + m'l) , (10) 


2 . S3T instability and emergence of finite ampli¬ 
tude large-scale structure 

The homogeneous equilibrium with no mean flow 

= = (13) 

is a fixed point of the S3T system (9) and (12) in the 
absence of hyperdiffusion (cf. Appendix B). The linear 
stability of the homogeneous equilibrium can be addressed 
by performing an eigenanalysis of the S3T system linearized 
about this equilibrium. The eigenfunctions in this case have 
the plane wave form 

5Z = , 5C = Cnm{x, 

(U) 

where x = Xi - X 2 , x = (a;i -I- X2)/2, y = y^ - y 2 , y = 
iUi + y 2 )l^, n and m are the x and y wavenumbers of the 
eigenfunction and tr = tr^ + is the eigenvalue with = 
Re(CT), ai = Im(tT) being the growth rate and frequency 
of the mode respectively. The eigenvalue a satisfies the 
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non-dimensional equation: 


~ pOO pOO 


{fak — fil) 

nfa{k'^ — 1 /.) + (fh? — TT?)k+l+ 

iP (kK^ 

-(k + n)k'^^ +{a + 2)k^k^ 


= {a + 1)7V^ — ifiP, (15) 

where L/ is a characteristic length scale, a = ajr and 
{n,m) = Lf(n,m) are the non-dimensional eigenvalue and 
wavenumbers respectively, e = e/(r^Ly) is the non-dimensional 
energy injection rate of the forcing, /3 = /3Ly/r is the non- 
dimensional planetary vorticity gradient, 

-| poo poo 

E{k,l) = —J J S(i,y)e-*“Midy , (16) 

is the Fourier transform of the forcing covariance, = 
k'^+P, Kg = {Ji+fiY + (J-\-fhY, iV^ = fc+ = k-\-n/2 

and = I + rh/2 (cf. Appendix B). For a forcing with the 
mirror symmetry 5(/c, —1) = 5(fc, 1) in wavenumber space 
and for h 0, the eigenvalues satisfy the relations: 

^(—n,m) and — ^{h,rh)i (1^) 

implying that the growth rates depend on |h| and \ fh\. As 
a result, the plane wave 5Z = cos{nx + my) and an array 
of localized vortices SZ = cos{nx) cos(mj/), have the same 
growth rate, despite their different structure. For zonally 
symmetric perturbations with n = 0, only the second rela¬ 
tion in (17) holds and (15) reduces to the eigenvalue relation 
derived by Srinivasan and Young (2012) for the emergence 
of jets in a barotropic /3-plane. 

We consider the case of a ring forcing that injects energy 
at rate e at the total wavenumber Kfi 

S(fc, 1) = 2Kf6{Vk^ + P- Kf), (18) 

and obtain the eigenvalues a by numerically solving (15). 
For small values of the energy input rate, < 0 for all 
wavenumbers and the homogeneous equilibrium is stable. 

At a critical ic the homogeneous flow becomes S3T unsta¬ 
ble and exponentially growing coherent structures emerge. 
The critical value, Sc, is calculated by first determining the 
energy input rate et(h, fa) that renders wavenumbers (n, fa) 
neutral (fTr(fi,m) = O) j and then by finding the minimum 
energy input rate over all wavenumbers: Sc = min(fi 
The critical energy input rate £c as a function of /3 is shown 
in figure 1. In addition, the corresponding critical zonos- 
trophy parameter Rjs = 0.7(ec,d^)^^^° which was used in 
previous studies to characterize the emergence and struc¬ 
ture of zonal jets in planetary turbulence (Galperin et al. 
2010), is shown as a function of /3 in figure 2. The absolute 


minimum energy input rate required is £c = 67 and occurs 
at Prriin = 3.5, while the minimum zonostrophy parameter 
required for the emergence of coherent flows is = 0.82 
and occurs for /3 — >■ 0. For (3 < I3min, the structures that 
first become marginally stable are zonal jets (with n = 0). 
The critical input rate increases as £c ~ for /3 —0 and 
the homogeneous equilibrium is structurally stable for all 
excitation amplitudes when f3 = 0. However, the structural 
stability for /3 = 0 is an artifact of the assumed isotropy of 
the excitation and the assumption of a barotropic flow. In 
the presence of even the slightest anisotropy (Bakas and 
loannou 2011, 2013b), or in the case of a stratified flow 
(Parker and Krommes 2015), zonal jets are S3T unstable 
and are expected to emerge even in the absence of /3. For 
j3 > Pminj the marginally stable structures are non-zonal 
and ic grows as ic /3^/^ for /3 —>■ oo. Since the critical 
forcing for the emergence of zonal jets (also shown in fig¬ 
ure 1), increases as £c ~ /3^ for /3 —>■ oo, for large values of 
j3 non-zonal structures first emerge and only at significantly 
higher e zonal jets are expected to appear. Investigation of 
these results with other forcing distributions revealed that 
the results for /3 ^ 1 are independent of the structure of 
the forcing (Bakas et al. 2015). 

The parameter regime of S3T instability is now related 
to the results of previous studies and to geophysical flows. 
Previous studies have identified a parameter regime which 
is distinguished by robust, slowly varying zonal jets as 
well as propagating, non-dispersive, non-zonal coherent 
structures (Galperin et al. 2010). This regime that is termed 
as zonostrophic, is in a region in parameter space in which 
the zonostrophy parameter is large {Rp > 2.5) and the 
scale kp = 0.5(/3^/e)^/® in which anisotropization of the 
turbulent spectrum occurs is sufficiently larger than the 
forcing scale {kp/Kf < 1/4). This regime is shown in 
figure 2 to be highly supercritical for all $. In addition, 
Bakas and loannou (2014) calculated indicative order of 
magnitude values of /3 and e for the Earth’s atmosphere 
and ocean as well as for the Jovian atmosphere. From these 
values we calculated the relevant zonostrophy parameter 
Rp and indicated the three geophysical flows in figure 2. 
We can see that all three cases are supercritical: the Jovian 
atmosphere is highly supercritical and is well within the 
zonostrophic regime, while the Earth’s atmosphere and 
ocean are slightly supercritical (at least within the context 
of the simplified barotropic model). 

We now examine the growth rate and dispersion prop¬ 
erties of the unstable modes for £ > £c and consider first 
the case /3 = 1, with e = 2£c. The growth rate of the maxi¬ 
mally growing eigenvalue, (7^, and its associated frequency 
of the mode, (7^, are plotted in figure 3(a) as a function of 
|h| and \fa\. We observe that the region in wavenumber 
space defined roughly by 0 < |fi| < 1/2, and 1/2 < |77i| < 1 
is unstable, with the maximum growth rate occurring for 
zonal structures (h = 0) with |7 ti| ~ 0.8. The frequency 
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Figure 1: The critical energy input rate Sc for structural 
instability (thick solid line) and the critical energy input rate 
for structural instability of zonal jets (solid line) as a function 
of The behavior of these critical values for 1 and 

/3 1 is indicated with the dashed asymptotes. In the light 

gray region only non-zonal coherent structures emerge, while 
in the dark gray region both zonal jets and non-zonal coherent 
structures emerge. The thin dotted vertical line (3 = ^rnin 
separates the unstable region: for (3 < Pmin zonal structures 
grow the most, whereas for /3 > ^rnin non-zonal structures 
grow the most. 



Figure 2: The critical zonostrophy parameter = 

0.7(£c/ 3^)^/^® for structural instability (thick line) and the 
corresponding critical parameter for structural instability of 
zonal jets (thin line) as a function of f3. The shaded region de¬ 
notes the zonostrophic regime for which both the inequalities 
> 2.5 and k^jKf < 1/4 are satisfied. The stars denote 
the position of the Earth’s atmosphere and ocean as well as 
the Jovian atmosphere in the {3 parameter space. 


Figure 3: Dispersion relation of the unstable modes for 
0 = 1 (panel a) and /3 = 10 (panel b). The contours show the 
growth rate dr and the shading shows the frequency dj of the 
unstable modes. For 0 0(1), stationary zonal jets are more 
unstable and for 0^1, westward propagating non-zonal 
structures are more unstable. For both panels, the energy 
input rate is e = 2dc. 


of the unstable modes is zero for zonal jet perturbations 
(n = 0) and non-negative for all other wavenumbers (n 0). 
Using the symmetries (17), this implies that real unstable 
mean flow perturbations SZ propagate in the retrograde 
direction if h ^ 0 and are stationary when n = 0. As e 
increases the instability region expands and roughly covers 
the sector 1/2 < |7V| < 1, with zonal structures having 
a larger growth rate compared to non-zonal structures, a 
result that holds for any e when P < Prnin- 

For P > Pmin the non-zonal structures have always 
larger growth rate. This is illustrated in figure 3(b), showing 
the growth rates and frequencies of the unstable modes 
for p = 10. For larger p values there is a tendency for 
the frequency of the unstable modes to conform to the 
corresponding Rossby wave frequency 


(^R 


Ph 

ffi -\- fh? ’ 


(19) 


a tendency that does not occur for smaller /3. A comparison 
between the frequency of the unstable modes and the Rossby 
wave frequency is shown in figure 4 in a plot of dildu. For 
slightly supercritical e, the ratio is close to one and the 
unstable modes satisfy the Rossby wave dispersion relation. 
At higher supercriticalities though, di departs from the 
Rossby wave frequency (by as much as 40% for the case of 
£ = SOEc shown in figure 4(b)). 
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Figure 4: Ratio of the frequency of the unstable modes 
ai over the corresponding frequency of a Rossby wave with 
the same wavenumbers at (a) e = 2ec and (b) e = 50ec 
when /3 = 100. Values of one denote an exact match with the 
Rossby wave frequency. 

3. Analysis of the eddy-mean flow dynamics under¬ 
lying jet formation 

In this section, we investigate the eddy-mean flow dy¬ 
namics leading to jet formation. These dynamics should 
have the property of directly channeling energy from the 
turbulent motions to the coherent flow without the pres¬ 
ence of a turbulent cascade. Previous studies have identified 
such mechanisms for the maintenance of zonal jets. Huang 
and Robinson (1998) showed that shear straining of the 
turbulent field by the jet produced upgradient momentum 
fluxes that maintained the jet against dissipation. A sim¬ 
ple case that clearly illustrates the physical picture for the 
mechanism of shear straining is to consider the evolution of 
eddies in a planar, inviscid constant shear flow. The eddies 
are sheared by the mean flow into thinner elliptical shapes, 
while their vorticity is conserved. For an elongated eddy 
this implies that the eddy velocities decrease and the eddy 
energy is transferred to the mean flow through upgradient 
momentum fluxes. This mechanism can operate when the 
time required for the eddies to shear over is much shorter 
than the dissipation time scale. The reason is that in this 
limit even the eddies with streamfunctions leaning against 
the shear that initially widen significantly gaining momen¬ 
tum, have the necessary time to shear over, elongate and 
surrender their momentum to the mean flow. Given that 
for an emerging jet the characteristic shear time scale is 
necessarily infinitely longer than the dissipation time scale, 
it needs to be shown that shear straining can produce up¬ 
gradient momentum fluxes in this case as well. In addition, 
previous studies have shown that shearing of isotropic eddies 
on an infinite domain does not produce any net momentum 
fluxes (Shepherd 1985; Farrell 1987; Holloway 2010) and 
should have no effect on the S3T instability (Srinivasan 


and Young 2012). Therefore another mechanism should be 
responsible for producing the upgradient fluxes in the case 
of an isotropic forcing. 

In order to investigate the eddy-mean flow dynamics 
underlying the S3T instability, we calculate the vorticity flux 
divergence that is induced when the statistical equilibrium 
(13) is perturbed by an infinitesimal coherent structure SZ. 
For an S3T unstable structure, the induced flux divergence 
tends to enhance the coherent structure SZ producing the 
positive feedback required for instability. So the goal of this 
section is to illuminate the eddy-mean flow dynamics leading 
to this positive feedback and to understand qualitatively 
why the homogeneous equilibrium is more stable for small 
and large values of (3. 

For zonal mean flows (9), (12) are simplified to: 
dtU = -dy {u'v')-rU = dy (20) 

and 

dtC=iA,+A2)C + E, ( 21 ) 

where 

A^ = -UA, - (/? - d^^yV)A-^d,, - r, (22) 

respectively. As a result the zonal mean flow is driven by the 
momentum flux divergence of the eddies. The perturbation 
in vorticity covariance SC that is induced by the mean 
flow perturbation SU can be estimated immediately by 
assuming that the system (20)-(21) is very close to the 
stability boundary, so that the growth rate is small. In this 
case the mean flow evolves slow enough that it remains in 
equilibrium with the eddy covariance, that is dSC/dt ~ 0. 
Bakas and loannou (2013b) showed that the ensemble mean 
momentum flux induced by an infinitesimal sinusoidal mean 
flow perturbation SU = esm(my), where e <C 1 (i.e the 
eigenfunction of (B4)), is equal in this case to the integral 
over time and over all zonal wavenumbers of the responses 
to all point excitations in the y direction: 

-I pOO pOO pOO 

S {u'v') = — / / u'v'(t)dtd^dk, (23) 

J—oo J—oo Jo 

where u'v'{t) is the momentum flux at time t produced by: 

G{k, y-i) = B{k)h{y - ^)Ukx+Uo{v-i)^ (24) 

The Green’s function G has the form of a wavepacket with 
an amplitude B(k) and a carrier wave with wavenumbers 
(fc, lo) that is modulated in the y direction by the wavepacket 
envelope h{y). The characteristics of the amplitude, the 
wavenumber and the envelope depend on the forcing char¬ 
acteristics, but in any case the calculation of the ensemble 
mean momentum fluxes is reduced to calculating the mo¬ 
mentum fluxes over the life cycle of wavepackets that are 
initially at different latitudes and then adding their relative 
contributions. 
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As the wavepacket propagates in the latitudinal direc¬ 
tion, its meridional wavenumber and frequency are going 
to change due to shearing by the mean flow and due to 
the change of the mean vorticity gradient /3 — Uyy. The 
resulting time variable momentum flux u'v’{t) can be calcu¬ 
lated using ray tracing. According to standard ray tracing 
arguments, the wave action is conserved along a ray (in the 
absence of dissipation) leading to the momentum flux: 

^{t) = -\B\'^AM{t)e-^^%{y - rim\ (25) 

where is the momentum flux of the 

carrier wave that determines the amplitude of the fluxes of 
the wavepacket and /*, ri{t) are the time dependent merid¬ 
ional wavenumber and position of the wavepacket respec¬ 
tively (Andrews et al. 1987). Because of the small amplitude 
of the mean flow perturbation SU, the wavenumber and 
position of the packet vary slowly on a time scale 0{et) 
compared to the dissipation time scale 1 /r and the domi¬ 
nant contribution to the time integral in (23) comes from 
small times. We can therefore seek asymptotic solutions of 
the form 


and slightly changes the amplitude of the fluxes as well as 
slightly speeds up or slows down the wavepacket. The sum 
of these two effects will produce the induced momentum 
fluxes. 

a. The limit of small scale wavepackets with a short propagation 

range 

In order to clearly illustrate the behavior of the eddy 
fluxes, we consider the limit of ,0 = j3Lf /r <C 1 , where Lf 
is the scale of the wavepackets and in addition we assume 
that the scale of the mean flow, 1 /m, is much larger than 
the scale of the wavepackets mLf <C 1. In this limit, the 
wavepackets are dissipated before propagating far from the 
initial position and the effect of the change in the mean 
vorticity gradient is higher order. As a result, Bakas and 
loannou (2013b) show that li and 771 decrease monotonically 
with time with rates independent of 5Uyy and proportional 
to the shear 5Uy{^) at the initial position 

h = -SUyiOkt , 771 = -mviO (28) 


— ^0 + + ■ ■ ■ I v{t) — C + cgt -l- £771 (t) -!-••• , (26) 


where cq = 2l3klo/{k'^ + 1^)^ is the group velocity in the 
absence of a mean flow and calculate the integral of u'v'{t) 
over time from the leading order terms. Substituting (26) 
in (25) we obtain: 


u'v'{t) = -\B\^AM{Q)e-^^%{y - ? - cot)|^ - 


-e\B\ 


dA 


M 


dlt 


li{t)e ^’'*|h(77-C-cot)|" 


- e\B\^AM{0)r,iit)e-^^^^\h{y - ^ - cot)\^ . 

ay 


u'v'i3 


(27) 


The first term, u'v'n, arises from the momentum flux pro¬ 
duced by a wavepacket in the absence of a mean flow. 
Because Am{ 0) = klo/{k‘^ + Iq)^ is odd with respect to 
wavenumbers, this term does not contribute to the en¬ 
semble averaged momentum flux when integrated over all 
wavenumbers and will be hereafter ignored. The second 
term, u'v's, arises from the small change in the amplitude 
of the flux Am over a dissipation time scale. The third 
term, u'v'p, arises from the small change in the position of 
the packet 77 compared to a propagating packet in the ab¬ 
sence of a mean flow. To summarize, the infinitesimal mean 
flow refracts the wavepacket due to shearing by the mean 
flow and due to the change of the mean vorticity gradient 


That is, the amplitude of the flux Am and the group velocity 
of the packets change only due to the shearing of the phase 
lines of the carrier wave according to the local shear. 

Consider in this limit the first term, u'v's, arising from 
the small amplitude change. Since the wavepacket is dis¬ 
sipated before it propagates away, we can ignore to first 
order propagation: 

^S=- e\B\^ h{t)e-^^^\h(y - ^ - c^t)\^ 

^\B\HUy{0kt(^^^ (29) 

so that the packet grows/decays in situ. Since the wave 
packet is rapidly dissipated, the integrated momentum flux 
over its life time will be given to a good approximation by 
the instantaneous change in the flux^ that is proportional 
to {dAM/dlt)ig. Figure 5 illustrates the amplitude of the 
momentum flux as a function of the angle 9t = arctan(lt/fc) 
of the phase lines of the carrier wave of the packet with the 
7 /-axis. It is shown that the momentum flux of a wavepacket 
with 1001 < 7 r /6 (that is with phase lines close to the merid¬ 
ional direction) excited in regions II or III, will increase 
within the dissipation time scale. Compared to an un¬ 
sheared wavepacket, this process leads to upgradient mo¬ 
mentum flux. The opposite occurs for waves excited in 
regions I and IV (with |0o| > 7 J'/ 6 ) that produce downgradi- 
ent flux, as their momentum flux decreases. 

We now consider the second term, u'v'js arising from 
the effect of propagation on the momentum flux. The group 

^occurring over the dissipation time scale 1/r that is incremental 
in shear time units 
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Figure 5: Amplitude of the momentum fluxes, Ajvf(t), of 
wavepackets as a function of the angle 6t = arctan(Zt/A:) 
between the phase lines of the central wave and the y-axis. 
The vertical lines separate the regions with \0t \ < tt/Q (II and 
III) and |i 9£| > 7r/6 (I and IV). 


velocity is given by Cg = 215 Am in this case and as a result a 
wavepacket starting in region III, will propagate towards the 
north (c.f. figure 5). Because shearing slows down the waves 
in region III [rji ~ —(dAM/dlt)), the wavepacket will flux 
its momentum from southern latitudes compared to when 
it moved in the absence of the shear flow. This is shown in 
figure 6(a) illustrating the distribution of momentum flux of 
an unsheared and a sheared perturbation whose amplitudes 
are constant. Figure 6(b) plots this difference, u'v'p, and 
shows that the flux is downgradient in this case. The same 
happens for waves excited in region II, while the waves 
excited in regions I and IV produce upgradient flux. 

The net momentum fluxes produced by an ensemble of 
wavepackets, will therefore depend on the spectral charac¬ 
teristics of the forcing that determine the regions (I-IV), in 
which the forcing has significant power. Bakas and loannou 
(2013b) show that for the isotropic forcing (18): 


<5 {u'v') = 

- J_ P 

“^ 7-00 


u'v' sd^dk 


3ej3‘^r d?5U 
‘52'kK'j dy^ 


1 

27r 



u'v'pd^dk 


(30) 


The first integral is zero, because the gain in momentum 
occurring for |6*o| < 7r/6 (waves excited in regions II, III) is 
fully compensated by the loss in momentum for |0o| > 7r/6 
(waves excited in regions I, IV) since for the isotropic forc¬ 
ing all possible wave orientations are equally excited. The 
net momentum fluxes are therefore produced by the u'v'p 
term and are upgradient, because the loss in momentum 
occurring for |0o| < 7r/6, is over compensated by the gain in 
momentum for |0o| > 7r/6. The momentum fluxes are also 


Figure 6: (a) Comparison of the momentum fluxes of an 
unsheared wavepacket excited in regions II (thick solid line) 
and III (solid line) to the momentum fluxes of a sheared 
wavepacket shown by the corresponding dashed lines, when 
only the change in propagation is taken into account. A 
snapshot of the fluxes at t = 0.2/r is shown. The planetary 
vorticity gradient is /3 = 0.1, the wavepacket has initial vortic- 

ity h{y) = e~^ , + Iq = 1, l^ol = tt/IO and |i?| = 1. (b) 

The difference in momentum fluxes between a sheared and an 
unsheared wavepacket calculated over their life cycle. 


proportional to the third derivative of SU yielding a hyper- 
diffusive momentum flux divergence that tends to reenforce 
the mean flow and is therefore destabilizing. These destabi¬ 
lizing fluxes are proportional to /3^ and as a result, the en¬ 
ergy input rate required to form zonal jets increases as l//3^ 
in this limit. It is worth noting that the first term integrates 
to zero only for the special case of the isotropic forcing, as 
even the slightest anisotropy yields a non-zero contribution 
from u'v's- For example consider the forcing covariance 
E!(a:i, a; 2 , j/i, j/ 2 ) = cos {k{xi — 0 : 2 )) that mim¬ 

ics the forcing of the barotropic flow by the most unstable 
baroclinic wave, which has zero meridional wavenumber. In 
this case the forcing that is centered at Zq = 0 in wavenum¬ 
ber space, injects significant power in a band of waves in 
regions II and III and therefore u'v's yields upgradient 
fluxes. 

b. The effects of the change in the mean vorticity gradient and 

the finite propagation range 

In order to take into account the effect of the change in 
the vorticity gradient, we retain higher order terms with 
respect to rnLf <C 1 in Zi and rji. In this case it can be 
shown that h decreases with time at a rate proportional 
to Uy{^) + Uyyy{^) (Balcas and loannou 2013b). Since the 
local shear and the local change in the vorticity gradient 
have different signs, the wavepacket is ’sheared less’ and as 
a result we expect reduced momentum fluxes compared to 
the limit discussed in section 4.2. Indeed, for the isotropic 













forcing: 


(d^5U 1 d^Su\ 

~32TTKj y dy^ ~ 4 ^ dy^ j ' 


(31) 


That is, the change in the mean vorticity gradient has a 
stabilizing effect. 

We finally relax the assumption that /3 <C 1. In this 
case, 1 1 and rji are affected by an integral shear and mean 
vorticity gradient over the region of propagation. For larger 
/3, the wavepacket will encounter regions of both positive 
and negative shear and as a result, the momentum fluxes 
that are qualitatively proportional to the integral shear over 
the propagation region will be reduced. In the limit /3 ^ I, 
the region of propagation is the whole sinusoidal flow with 
consecutive regions of positive and negative shear and the 
integral shear along with the fluxes will asymptotically tend 
to zero. As a result, the energy input rate required for 
structural instability of zonal jets increases with (3 in this 
limit. 


4. Equilibration of the S3T instabilities 

We now investigate the equilibration of the instabilities 
by studying the S3T system (9), (12) discretized in a doubly 
periodic channel of size 27r x 27r. We approximate the 
monochromatic forcing (18), by considering the narrow 
band forcing 


s(fc,o 


r 1, for A/l < AKf 

AKf \ 0, for \Vk‘^ + P-Kf \ > AKf ’ ^ ’ 


where fc, I assume integer values, that injects energy at rate 
e in a narrow ring in wavenumber space with radius Kf and 
width AKf. We consider the set of parameter values (3 = 10, 
r = 0.01, v = 1.19 • 10“®, Kf = 10 and AKf = 1, for which 
P = 100. The integration is therefore in the parameter 
region of hgure 1 in which the non-zonal structures are 
more unstable than the zonal jets. The growth rates of the 
coherent structures for integer values of the wavenumbers, 
n and m are calculated from the discrete version of equation 
(15) obtained by substituting the integrals with sums over 
integer values of the wavenumbers (Bakas and loannou 
2013a). 

We first consider the supercritical energy input rate 
e = dCc- For these parameters only non-zonal modes are un¬ 
stable, with the perturbation with {n,m) = (1,5) growing 
the most. At t = 0, we introduce a small random perturba¬ 
tion, whose streamfunction is shown in hgure 7(a). After 
a few e-folding times, a harmonic structure of the form 
Z = cos(a:) cos(5?/) dominates the large-scale how. The 
energy of this large scale structure shown in hgure 7(b), 
increases rapidly and eventually saturates. At this point 
the large-scale how gets attracted to a traveling wave hnite 
amplitude equilibrium structure (cf. hgure 7(c)) close in 



Figure 7: Equilibration of the S3T instabilities, (a) Stream- 
function of the initial perturbation, (b) Energy evolution 
of the initial perturbation shown in panel (a) as obtained 
from the integration of the S3T equations (9) and (12) 
(dashed line) and from the integration of the ensemble quasi- 
linear (EQL) system (4)- (3) with A^ens = 10 (solid line) and 
IVens = 100 (dash-dotted line) ensemble members that is dis¬ 
cussed in section 6. (c) Snapshot of the streamfunction '^eq 
of the traveling wave structure and (d) Hovmoller diagram 
of '^eq{x, y = 7r/4, t) for the finite equilibrated traveling wave. 
The thick dashed line shows the phase speed obtained from 
the stability equation (15). The energy input rate is i = 4ec 
and jS = 100. 


form to the harmonic Z = cos(x) cos(5?/) that propagates 
westward. This is illustrated in the Hovmoller diagram of 
'0(a;, y = 7r/4, t) shown in 7(d). The sloping dashed line in 
the diagram corresponds to the phase speed of the traveling 
wave, which is found to be approximately the phase speed 
of the unstable {n,m) = (1,5) eigenmode. 

Consider now the energy input rate e = lOec- While the 
maximum growth rate still occurs for the (|n|, \m\) = (1, 5) 
non-zonal structure, zonal jet perturbations are unstable 
as well. If the S3T dynamics are restricted to account only 
for the interaction between zonal flows and turbulence by 
employing a zonal mean rather than an ensemble mean, 
an infinitesimal jet perturbation will grow and equilibrate 
at finite amplitude. To illustrate this we integrate the 
S3T dynamical system (20)-(21) restricted to zonal flow 
coherent structures. The energy of the small zonal jet 
perturbation 5Z = O.lcos(42/) is shown in figure 8 to grow 
and saturate at a constant value and the streamfunction 
of the equilibrated jet is shown in the left inset in figure 
8. However, in the context of the generalized S3T analysis 
that takes into account the dynamics of the interaction 
between coherent non-zonal structures and jets, we find 
that these S3T jet equilibria can be saddles: stable to zonal 
jet perturbations but unstable to non-zonal perturbations. 
To show this, we consider the evolution of the same jet 
perturbation 5Z = 0.1cos(4j/) under the generalized S3T 
dynamics (9), (12) and find that the flow follows the zonally 
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Figure 8: Energy evolution of an initial jet perturbation 
SZ = 0.1 cos(4j/) for the zonally restricted S3T dynamics 
(20)-(21) (dashed line) and the generalized S3T dynamics (9), 

(12) (thin line). The insets show a snapshot of the mean 
flow streamfunction at t = 700 (left) and the streamfunction 
of the equilibrated structure at i = 3500 (right) under the 
generalized S3T dynamics. The parameters are e = lOec and 
/3 = 100 . 

restricted S3T dynamics and equilibrates to the same finite 
amplitude zonal jet (cf. figure 8). At this point we insert a 
small random perturbation to the equilibrated flow. Soon 
after, non-zonal undulations grow and the flow transitions 
to the stable Z = cos(a;) cos(5y) traveling wave state that 
is also shown in figure 8. As a result, the finite equilibrium 
zonal jet structure is S3T unstable to coherent non-zonal 
perturbations and is not expected to appear in non-linear 
simulations despite the fact that the zero flow equilibrium 
is unstable to zonal jet perturbations. We will elaborate 
more on this issue in the next section. 

Finally, consider the case e = 30£c- At this energy input 
rate, the finite amplitude non-zonal traveling wave equilibria 
become S3T unstable. To show this, we consider the non- 
zonal traveling wave equilibrium obtained by the evolution 
of the small non-zonal perturbation 5Z = 0.01 cos(x) cos(5y) 
to the homogeneous state that is shown at the left inset 

in figure 9 and impose a small random zonal perturbation. 

_2 

The evolution of the zonal energy = (l/2)t/ , where the 
overbar denotes a zonal average, is shown in figure 9. After 
an initial transition period, the zonal perturbations grow 
exponentially and the flow transitions to the jet equilibrium 
state shown at the right inset in figure 9. Note however, 
that the jet equilibrium structure is not zonally symmetric. 
This is a new type of S3T equilibrium: it is a mix between 
a zonal jet and a non-zonal traveling wave with the same 
meridional scale. These mixed equilibria appear to be the 
attractors for larger energy input rates as well. This is 
illustrated in figure 10 showing the structure of the mixed 
equilibrium at £ = 50£c. The equilibrium structure consists 
of a large amplitude zonally symmetric jet with larger scale 


Figure 9: Zonal energy evolution of a random zonal pertur¬ 
bation imposed on the non-zonal traveling wave equilibrium 
shown in the left inset. The streamfunction of the equilibrated 
structure is shown in the right inset. The energy input rate 
is £ = SOsc and ^ = 100. 

compared to the mixed state in figure 9. Embedded in it 
are non-zonal vortices with the same meridional scale and 
with about 14% the energy of the zonal jet. These vortices 
that are shown in figure 10(b) to have approximately the 
compact support structure T = cos(x) cos(4y) propagate 
westward as shown in the Hovmoller diagram in figure 10(c). 

5. Comparison to ensemble mean quasi-linear and 
non-linear simulations 

a. Comparison to an ensemble of quasi-linear simulations 

Within the context of the second order cumulant closure, 
the S3T formulation allows the identification of statistical 
turbulent equilibria in the infinite ensemble limit, in which 
the fluctuations induced by the stochastic forcing are av¬ 
eraged to zero. However, these S3T equilibria and their 
stability properties are manifest even in single realizations 
of the turbulent system. For example, previous studies us¬ 
ing S3T obtained zonal jet equilibria in barotropic, shallow 
water and baroclinic flows in close correspondence with 
observed jets in planetary flows (Farrell and loannou 2007, 
2008, 2009b, a). In addition, previous studies of S3T dynam¬ 
ics restricted to the interaction between zonal flows and 
turbulence in a /3-plane channel showed that when the en¬ 
ergy input rate is such that the zero mean flow equilibrium 
is unstable, zonal jets also appear in the non-linear simu¬ 
lations with the structure (scale and amplitude) predicted 
by S3T (Srinivasan and Young 2012; Constantinou et al. 
2014). 

A very useful intermediate model that retains the wave- 
mean flow dynamics of the S3T system while relaxing the 
infinite ensemble approximation is the quasi-linear system 
(5)-(6). Under the ergodic assumption, this can be inter- 


10 



































Figure 10: Mixed zonal jet-traveling wave S3T equilibrium 
for € = 50£c and /3 = 100. (a) Snapshot of the streamfunction 
'^eq of the equilibrium state, (b) Contour plot of the non-zonal 
component ’J'eg — '^eq of the equilibrium structure, where the 
overline denotes a zonal average, (c) Hovmoller diagram of 
^eqi^^y = 7r/4,f) for the equilibrated structure. 


preted as an ensemble of quasi-linear equations (EQL) in 
which the ensemble mean can be calculated from a finite 
number of ensemble members. Its integration is done as 
follows. A pseudo-spectral code with a 128 x 128 resolution 
and a fourth order Runge-Kutta scheme for time stepping 
is used to integrate (5)-(6) forward. At each time step, Nens 
separate integrations of (6) are performed with the eddies 
evolving according to the instantaneous flow. Then the 
ensemble average vorticity flux divergence is calculated as 
the average over the N^ns simulations and (5) is stepped 
forward in time according to those fluxes. The EQL system 
reaches a statistical equilibrium at time scales of the order 
of tf,q ~ 0(l/r) and the integration was carried on until 
t = lOOteq in order to collect accurate statistics. 

We choose the same parameter values as in the S3T 
integrations in section 5 (/3 = 10, r = 0.01, = 1.19 • 10“®, 

Kf = 10 and AKf = 1). For these parameters (/3 = 
100), S3T predicts the emergence of propagating non-zonal 
structures when the energy input rate exceeds the critical 
threshold ec^ and the emergence of mixed zonal jet-traveling 
wave states when the finite amplitude traveling wave states 
become structurally unstable to zonal jet perturbations. 
In order to examine whether the same bifurcations occur 
in the EQL system, we consider two indices that measure 
the power concentrated at scales larger than the scales 
forced. The first is the zonal mean flow index defined as in 
Srinivasan and Young (2012), as the ratio of the energy of 
zonal jets with scales larger than the scale of the forcing 
over the total energy 


zmf = 


J2l:l<Kf-AKf 0 

Ekimi) 


(33) 


Figure 11: The zmf and nzmf indices defined in (33) and 
(35) respectively, as a function of energy input rate e/sc and 
the zonostrophy parameter for the non-linear (NL) inte¬ 
grations and an ensemble of quasi-linear (EQL) integrations 
(dashed line) with Nens = 10 ensemble members as described 
in section 6. The critical value Sc = 8.4 • 10~® is the energy 
input rate at which the S3T predicts structural instability 
of the homogeneous turbulent state. Zonal jets emerge for 
with £ni = IbSc- The parameters are /3 = 10, 
r = 0.01, u = 1.19-10~® and the forcing is an isotropic ring in 
wavenumber space with radius Kf = 10 and width AKf = 1. 


where 


E{kJ) 



icr \, ) 


dt 


(34) 


is the time averaged total energy power spectrum of the 
flow at wavenumbers {k,l). The second is the non-zonal 
mean flow index defined as the ratio of the energy of the 
non-zonal modes with scales larger than the scale of the 
forcing over the total energy: 

^ J2kl:K<Kf-AKf 

EuEikJ) 

If the structures that emerge are coherent, then these indices 
quantify their amplitude. Figure 11 shows both indices as 
a function of the energy input rate e and as a function 
of the corresponding values of the zonostrophy index Rjs 
for EQL simulations with Nens = 10 members. The rapid 
increase of the nzmf index for e > Sc (corresponding to 
RjS > 1.55), illustrates that this regime transition in the flow 
predicted by S3T with the emergence of non-zonal structures 
manifests in the quasi-linear dynamics as well. We now 
consider the case e = isc in detail in which the traveling 
wave structure Z = cos(a:) cos(5?/) is maintained in the S3T 
integrations. We observe, that the S3T equilibria manifest 
in the EQL simulations with the addition of some ’thermal 
noise’ due to the stochasticity of the forcing that is retained 
in this system. This is illustrated in figure 7b showing the 
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Figure 12: Snapshot of the mean streamfunction 'F at sta¬ 
tistical equilibrium obtained from the ensemble mean quasi- 
linear simulations with Nena = 10 members for e = Asc (panel 
a) and e = SOec (panel b). The parameters are as in figure 11. 


energy growth of the coherent structure for N^ns = 10 and 
IVens = 100. The energy of the coherent structure in the 
EQL integrations fluctuates around the values predicted 
by the S3T system with the fluctuations decreasing as 
1 /V^ens- However, even with only 10 ensemble members we 
get an estimate that is very close to the theoretical estimate 
of the infinite ensemble members obtained from the S3T 
integration. The structure of the traveling wave equilibrium 
in the quasi-linear simulations shown in figure 12(a) and its 
phase speed (not shown) are also in very good agreement 
with the corresponding structure and phase speed obtained 
from the S3T integration. 

The second transition in which zonal jets emerge is 
more intriguing. While the homogeneous equilibrium is 
structurally unstable to zonal jets when Ssz = 5.2ec, the 
finite amplitude zonal jet equilibria are structurally unsta¬ 
ble and the flow stays on the attractor of the non-zonal 
traveling wave equilibria (cf. figure 8). When e > the 
non-zonal traveling wave equilibria become S3T unstable 
while at these parameter values the S3T system has mixed 
zonal jet-traveling wave equilibria which are stable (cf. fig¬ 
ure 10). The rapid increase in the zmf index with the 
concomitant rapid decrease in the nzmf index shown in 
figure 11, illustrates that this regime transition manifests in 
the EQL system as well with similar mixed zonal-traveling 
wave states appearing. The structure of the mixed zonal 
jet-traveling wave equilibrium for e = bOSc is shown in 
12(b) and similar to the S3T equilibrium in hgure 10, it 
consists mainly of 4 zonal jets and the compact support 
vortices Z ~ cos(a:) cos(4?/) embedded in the jets. We there¬ 
fore conclude that the EQL system accurately captures the 
characteristics of the emerging structures. 

b. Comparison to non-linear simulations 

In order to compare the predictions of S3T to the non¬ 
linear simulations, we solve (1) with the narrow band forcing 
(32) on a doubly periodic channel of size 27r x 27r using the 
same pseudospectral code as in the EQL simulations and 
the same parameter values. Figure 11 shows the nzmf and 
zmf indices as a function of the energy input rate e for the 



Figure 13: Time averaged energy power spectra, log(.B(fc, i)), 
obtained from the non-linear (NL) simulation of (1) at e/sc = 
4 (panel a) and e/ec = 50 (panel c). Hovmoller diagram of 
tp{x,y = nl4,t) at e/ec = 4 (panel b) and e/ec = 50 (panel 
d). The thick dashed lines correspond to the phase speed 
obtained from the eigenvalue relation (15). 


NL simulations. The rapid increase in the nzmf index for 
e > £c shows that the non-linear dynamics share the same 
bifurcation structure as the S3T statistical dynamics. In 
addition, the stable S3T equilibria are in principle viable 
repositories of energy in the turbulent flow and the non¬ 
linear system is expected to visit their attractors for finite 
time intervals. Indeed for e = the pronounced peak 
at (|fc|, \l\) = (1,5) of the time averaged power spectrum 
shown in figure 13(a) illustrates that the traveling wave 
equilibrium with (|fc|, |(|) = (1,5) that emerges in the S3T 
integrations, is the dominant structure in the NL simula¬ 
tions. Comparison of the energy spectra obtained from 
the EQL and the NL simulations (not shown), reveals that 
the amplitude of this structure in the quasi-linear and in 
the non-linear dynamics almost matches. Remarkably, the 
phase speed of the S3T traveling wave matches with the 
corresponding phase speed of the (|A:|, |(|) = (1,5) struc¬ 
ture observed in the NL simulations, as can be seen in the 
Hovmoller diagram in figure 13(b). Such an agreement in 
the characteristics of the emerging structures between the 
EQL and NL simulations occurs for a wide range of energy 
input rates as can be seen by comparing the nzmf indices in 
figure 11. As a result, S3T predicts the dominant non-zonal 
propagating structures in the non-linear simulations, as well 
as their amplitude and phase speed. 

We now focus on the second regime transition with the 
emergence of zonal jets. The increase in the zmf index in 
the NL simulations for e > Sni that is shown in figure 11, 
indicates the emergence of jets roughly at the bifurcation 
point of the S3T and EQL simulations. However, the energy 
input rate threshold for the emergence of jets is larger in the 
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NL simulations compared to the corresponding EQL thresh¬ 
old. This discrepancy possibly occurs due to the fact that 
the exchange of instabilities between the mixed jet-traveling 
wave equilibria and the pure traveling wave equilibria de¬ 
pends on the equilibrium structure \Z^ ^ C^]. Small changes 
for example in that might be caused by the eddy-eddy 
terms neglected in S3T can cause the exchange of insta¬ 
bilities to occur at slightly different energy input rates. It 
was shown in a recent study that when the effect of the 
eddy-eddy terms is taken into account by obtaining 
directly from the nonlinear simulations, the S3T stability 
analysis performed on this corrected equilibrium states ac¬ 
curately predicts the energy input rate for the emergence of 
jets in the nonlinear simulations (Constantinou et al. 2014). 
The power spectrum obtained from the NL simulations for 
e = 50ec shows an energy peak at (fc, \l\) = (0,4) with sec¬ 
ondary power peaks at (|A:|, |/|) = (1,4) and (|/c|, \l\) = (1, 5)) 
(of approximately 12% of the energy in the zonal jet each). 
The Hovmoller diagram of the streamfunction shown in 
figure 13(d) reveals that the dominant non-zonal structures 
in the NL simulations propagate in the retrograde direction. 
As a result the mixed S3T equilibrium of figure 10 manifests 
in the NL simulations. Note however, that the phase speed 
calculated from the diagram is different than the phase 
speed of the (|fc|,|Z|) = (1,4) structure in figure 10. At 
larger energy input rates the zonal jets have typically larger 
scales due to jet merging and coexist with energetically sig¬ 
nificant westward propagating non-zonal structures having 
an energy between 10 — 50% of the jet energy and scales 
(|fc|,|Z|) = (1 ,to), where m is the number of jets in the 
channel. Such an agreement again holds for a wide range 
of energy input rates, as the zmf indices obtained from the 
EQL and the NL simulations indicate. In summary, S3T 
predicts the characteristics of both non-zonal propagating 
structures and of zonal jets in the non-linear simulations. 

c. Zonostrophic regime 

S3T and the corresponding ensemble quasi-linear sys¬ 
tem were obtained by ignoring the eddy-eddy non-linear 
interactions. Therefore the question arises on whether the 
predictions of S3T are useful in the zonostrophic regime. In 
this regime, which is highly supercritical with respect to S3T 
instability of the homogeneous equilibrium (cf. figure 2), 
maintenance of zonal jets and zonons were interpreted by 
previous studies to arise from an inverse energy cascade 
(Galperin et al. 2010), a highly non-linear process, which is 
absent in S3T. According to this interpretation, the turbu¬ 
lent energy cascades isotropically toward large scales until 
it reaches kp. At this scale the cascade becomes anisotropic 
and most of the energy is channelled into the zonal flows. 
To illustrate this, the time averaged energy power spec¬ 
tra E{k,l) are typically split between the zonal spectra 
Ezil) = E{k = 0,1) and the residual Eji{k,l) = E — E^. 
The zonal and residual spectra calculated from NL inte¬ 


grations in the zonostrophic regime (AT/ = 60, AK = 1, 
/3 = 42, r = 0.01, e = 0.0065) are shown in figure 14. Up 
to the scale kp, the residual spectra follow the Kolmogorov 
K~^l^ law in accordance with an isotropic cascade assump¬ 
tion. At this scale, the cascade is anisotropized and the 
residual spectra steepen. However, most of the energy is in 
zonal scales with the zonal spectra following a much steeper 
K-^ law. 

The residual and the zonal spectra obtained from an 
EQL simulation with A^ens = 10 for the same parameters, 
are also shown in figure 14. The residual spectra follow a 
shallower than slope for K > kp, while they steepen 

after kp and reach a lower peak with respect to the cor¬ 
responding spectra from the NL simulations. In addition, 
the residual part of the spectrum corresponds mainly to 
incoherent motions for scales with K > kp. This is revealed 
by taking into account only the spectra of the coherent 
part of the flow and calculating the residual spectrum that 
is also shown in figure 14. For most of the scales, it is 
at least one to two orders of magnitude lower than the 
corresponding residual spectrum when both coherent and 
incoherent motions are taken into account and only the 
non-zonal structures with large scales (close to the energy 
peak) appear to be coherent. Failure of the EQL simula¬ 
tions to exactly reproduce the slope of the incoherent 

turbulent motions is not surprising, since the inverse energy 
cascade that is absent in the EQL simulations is essential 
for this part of the spectrum. The energetically important 
part however, which contains the large-scale energetic waves 
is captured by S3T. The zonal spectra obtained from the 
EQL simulations follow the same K~^ law and peak at the 
same scale compared to the NL simulations but the peak 
has a larger amplitude. As is argued in section 5.2.2, the 
steep power law is an artifact of the shape of the strongly 
forced jet which is characterized by near discontinuity in 
the shear at the maxima of the prograde jets. 

So to summarize, the scale and the shape of the domi¬ 
nant jet structure, as well as the scale of the most energetic 
coherent non-zonal structures are accurately captured by 
the EQL simulations, while the eddy-eddy interactions ne¬ 
glected in the EQL simulations set the proper scaling for 
the tail of the spectrum that consists of incoherent turbu¬ 
lent motions and change the partition between the energy 
of the jet (that is overestimated in the EQL simulations) 
and that of the non-zonal large-scale structures (that is 
underestimated in the EQL simulations). 

6. Summary — Discussion 

This article addressed the emergence of coherent struc¬ 
tures in barotropic /3-plane turbulence using the tools of 
Stochastic Structural Stability Theory (S3T), a statistical 
theory that expresses the statistics of the turbulent flow 
dynamics as a systematic cumulant expansion truncated 
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Figure 14: Residual (blue) and zonal (red) energy spectra 
for the NL (solid) and EQL (dashed) simulations in the zonos- 
trophic regime. Also shown are the residual spectra from the 
EQL simulations when only the coherent motions are taken 
into account (blue dotted line). The parameters are Kf = 60, 

AK = 1, j3 = 42, r = 0.01, e = 0.0065, for which ^ = 70, 

£ = 7308ec, kp = 12.9 and Rp = 2.5. Lines (thin dashed) 
of slope K~^l^ and K~^ are also plotted for reference. The 
pseudo-spectral code was run with a 512 X 512 resolution 
and the exponential filter of Smith et al. (2002) instead of 
hyper-diffusion. 

at second order. With the interpretation of the ensemble 
average as a Reynolds average over the fast turbulent eddies 
adopted in this contribution, the second order cumulant 
expansion results in a closed, non-linear dynamical system 
that governs the joint evolution of slowly varying, spatially 
localized coherent structures with the second order statistics 
of the rapidly evolving turbulent eddies. The fixed points 
of this autonomous, deterministic non-linear system define 
statistical equilibria, the stability of which are amenable 
to the usual treatment of linear and non-linear stability 
analysis. 

The linear stability of the homogeneous S3T equilibrium 
with no mean velocity was examined analytically for the case 
of an isotropic random stirring that sustains the turbulence 
in the barotropic flow. Structural instability was found to 
occur for perturbations with smaller scale than the forcing, 
when the energy input rate e = eK^j/r^ is larger than a 

certain threshold ic that depends on ,5 = j3/{rKf). It was 
found that when /? is small or order one, the maximum 
growth rate occurs for stationary zonal structures, while for 
large /3, westward propagating non-zonal grow the most. 

The eddy-mean flow dynamics underlying the S3T in¬ 
stability of zonal jets was then studied in detail. It was 
shown that close to the structural stability boundary, the 
dynamics can be split into two competing processes. The 
first, which is shearing of the eddies by the local shear 
described by Orr dynamics in the /? plane, was shown to 
lead to jet forming upgradient momentum fluxes acting 


exactly as negative viscosity for an anisotropic forcing and 
as negative hyperviscosity for isotropic forcing. The second, 
which is momentum flux divergence resulting from lateral 
wave propagation on the nonuniform local mean vorticity 
gradient, was shown to lead to jet opposing downgradient 
fluxes acting as hyperdiffusion. 

The equilibration of the unstable, exponentially growing 
coherent structures for large /3 was then studied through 
numerical integrations of the S3T dynamical system. When 
the forcing amplitude is slightly supercritical, the finite am¬ 
plitude traveling wave equilibrium has a structure close to 
the corresponding unstable non-zonal perturbation with the 
same scale. When the forcing amplitude is highly supercrit¬ 
ical, the instabilities equilibrate to mixed states consisting 
of strong zonal jets with smaller amplitude traveling waves 
embedded in them. 

The predictions of S3T were then compared to the re¬ 
sults obtained from direct numerical simulations of the 
turbulent dynamics. The critical threshold above which 
coherent non-zonal structures are unstable according to the 
stability analysis of the S3T system was found to be in 
excellent agreement with the critical value above which non- 
zonal structures acquire significant power in the non-linear 
simulations. The scale, phase speed and amplitude of the 
dominant structures in the non-linear simulations were also 
found to correspond to the structures predicted by S3T. 
In addition, the threshold for the emergence of jets, which 
is identified in S3T as the energy input rate at which an 
S3T stable, finite amplitude zonal jet equilibrium exists, 
was found to roughly match the corresponding threshold 
for jet formation in the non-linear simulations, with the 
emerging jet scale and amplitude being accurately obtained 
using S3T. Such a good agreement between the predictions 
of S3T and direct numerical simulations, holds not only 
close to the bifurcation point for the emergence of coherent 
structures but also in the regime of zonostrophic turbulence. 
Consequently, these results provide a concrete example that 
large-scale structure in barotropic turbulence, whether it is 
zonal jets or non-zonal coherent structures, emerges and is 
sustained from systematic self-organization of the turbulent 
Reynolds stresses by spectrally non-local interactions and 
in the absence of a turbulent cascade. 

APPENDIX A 

Boundedness of the solutions and invariants of the 
S3T equations 

The S3T system in the absence of forcing and dissipation 
has similar quadratic invariants with the nonlinear system. 
Further, the solutions of the S3T equations remain bounded 
for all times. That is, the sum of the enstrophy of the ensem¬ 
ble mean over the domain, Hm = 1/2/ Z^dxdy, and the 
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eddy enstrophy over the domain, Hp = 1/2 J 
is conserved. Similarly, the sum of the energy of the en¬ 
semble mean, E^n = 1/2 / (?7^ -I- V^)dxdy, and the eddy 
energy, Ep = 1/2 J{A^^C)xi=x 2 dxdy, is also conserved. 
We show this by first multiplying (9) (in the absence of 
hyper-diffusion) by Z to obtain: 

dtVm + Udxr]m + VdyT]m + jdVZ = -ZV ■ (u'C') - 2r?7m, 

(Al) 

where rjm = Z'^ 12 is the enstrophy density of the ensem¬ 
ble mean. Integrating by parts and using the continuity 
equation we rewrite the advection terms as: 

Uda^ym + Vdy^m = d^iU ym) + dyiVyrn)- (A2) 

Writing Z = d^V — dyll and using again the continuity 
equation we have: 

1/2 + y2 

ZV = d ,—I- dyiUV), (A3) 

and (Al) becomes: 


In this Appendix we derive the dispersion relation (15), 
which determines the stability of zonal as well as non-zonal 
perturbations in homogeneous turbulence. We follow closely 
the treatment of Srinivasan and Young (2012). We first 
rewrite (9), (12) in terms of the variables x = Xi — X 2 , 
(l/2)(a:i-fa:2), y = yi-y 2 and y = (l/2)(yi-^^ 2 )- The 
derivatives transform into this new system of coordinates to 

= (l/2)^+(-l)*+la£, dy^ = il/2)dy + {-lY + ^dy, A, = 

A + (_l/4)A-f (-l)^+ia?^ + (-l)*+iaL, with A = dl, + dlp 
and A = + d^. It is also convenient to introduce the 

streamfunction covariance S{x,x,y,y) = {ijji'YY)^ which is 
related to C{x,x,y,y) via: 


c = (c;c 2 ) = (AiV''iA2 V'^) = A1A25 


A- 


A -r^ 


s, 


(Bl) 


where T = i9|j -|- d'^y. Equations (9), (12) then become in 
the absence of hyper-viscosity (^ = 0): 


dtym+^ ■{^ym)+l3dxem-lddy{UV) = -ZV-(u'C')-2r??m, 

(A4) 

where = (1/2)(U^ + V"^) is the energy density of the 
ensemble mean. Similarly it can be shown from (12), 
that the ensemble mean of the eddy enstrophy density 
yp = ( 1 / 2 )Cxi=x 2 evolves (in the absence of hyper-diffusion) 
according to: 


dt + Ud^+Udi + Vdy + Vdy 


c+ 


(/? + Zy)dx Zydx — Zxdy — ZxBy 


A + ) 5- 


2 (/? -|- Zy)dx + —Zydx ~ ‘2.Zxdy — -^Zxdy 


= -2rC' + S, 


rs' = 


(B2) 


dtyp + V • (Uryp) + (3 [^^^(ep) - dy {u'v')] + 

+ {u'C') dxZ + (v'C) dyZ = yf- 2ryp, (A5) 

where Sp = (l/ 2 )(A/^C')xi=x 2 is the ensemble mean of the 
eddy energy density and yj = (l/ 2 )Sxi=x 2 is the enstrophy 
density of the forcing. Adding (A4) and (A5) we obtain the 
equation for the evolution of the total enstrophy density 

»? = ??p + dm- 

{dt+2r)y-yf = -V ■{\Jy)-l3dx{epEem)+Y^dy{UV+{u'v'))- 

(A6) 

Integrating (A6) over the horizontal domain, the terms on 
the rhs of (A6) integrate to zero and the total enstrophy 
El = Hjn + Hp = f (yxn + yp)dxdy evolves according to: 


(a* + u • V) z + /3E = {dly - a?^)rAU=5=o - rz, (b3) 

where (i7,F) = (l/2)(17i + U 2 , W + F), (F, F) = (C/i - 

U 2 ,Vi-V 2 ), {Zx,Zy) = {l/2){dx2 +dx 2 ,dy^ + c>y 2 )Z and 

{Zx,Zy) = {dx -2 - dx 2 ,dy, - dyY)Z. 

The forcing covariance 5 is homogeneous and as a result 
it depends only on the difference coordinates, x and y. It 
can then be readily shown from (B2)-(B3), that the state 
with no coherent flow {U^ = = Z^ = 0) and with 

the homogeneous vorticity covariance C^{x,y) = S/(2r) 
(implying also that the streamfunction covariance is 
homogenous) is a fixed point of the S3T system. The 
stability of this homogeneous equilibrium, can be addressed 
by first linearizing the S3T system about it: 


dtH ^Hf- 2rH, {AT) 

where Hf is the total enstrophy imparted by the forcing. 
As a result, the enstrophy is bounded and has the value 
= Hf/{2r) at steady state. Similarly, it can be shown 
that the total energy E = Em + Ep is bounded. 


dtdC = -( 5Udx + 5Vd, 




-(5Z,dx-5ZxdAAS^- 


A + yA 

4 

dt5Z = -m + {dly - dUr6S\x=y=o - rdZ, 


dx-2Tdx>6S -2rSC, (B4) 


(B5) 


APPENDIX B 

Calculation of the dispersion relation and its 
properties 


where SZ, 5U, 5V, 5Zx^ 5Zy, 5C and 5S are small per¬ 
turbations in the ensemble mean vorticity, velocities and 
vorticity gradients and in the eddy vorticity and stream- 
function covariances respectively, and then performing an 
eigenanalysis of the linearized equations (B4)-(B5). 
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We consider a harmonic vorticity perturbation of the 
form 8Z = for which: 

[5UJVJZ,,6Zy] = 

r, r ^ 1 . /na; my\ ,r,^.irnT, rrt 

= -2 It + 2 j 

(B6) 

with 7V^ = n? + m?. Taking the same form for the stream- 

function covariance perturbation SS = Snm{x, 

and inserting it in (B4)-(B5) along with (B6) yields: 

2 


(cr -I- 2r) 


A-- A 


2il3A+dx — inp ( A — 




Sr}.7T}. — 


= — sm 
7V2 


^ ^ - ndy) (A -h N‘^)KS^ 


(B7) 

- (cr -h r)N‘^ + inl3 = {mds, - ndy) A+Snm\x=y=o, 

(B8) 

where A+ = ndx + mdy and C® = S/2r = A'^S^ is the 
equilibrium vorticity covariance with zero mean flow. 
Defining the Fourier transform of Snmix,y) by 

-I poo pOO 

Snm(.k,l) = — / 5„™(i,y)e-*“Midy, (B9) 

J — oo J — oo 


we obtain from (B7) that the Fourier component Snm sat¬ 
isfies: 

~ {mk--nl_)Kl{Kl/N‘^-l)S^{k_,l_) 

^ il3{k_Kl - k+Kl) + {a + 2r)KlKl 

{mk+-nl+)Kl{Kl/N‘^ - l)S’^{k+,l+) 
il3{k-Kl-k+Kl) + {a + 2r)KlKl ’ 

(BIO) 

with k± = k ± n/2, l± = I ± m/2, = k'^ + and 

= k"^ + P. = Ej(2rK‘^) is the Fourier transform of 
, and S is the Fourier transform of S. In addition, (B8) 
becomes: 


Equation (Bll) can be further simplified by noting that 
because the choice of xi and X 2 is arbitrary, the forcing co- 
variance satisfies the exchange symmetry S(a;i,^ 2 , J/i, 2 / 2 ) = 
S(x 2 , Xi,y 2 ,yi)- In terms of the new variables, the exchange 
symmetry is written as 'E.{x,x,y,y) = 'A{—x,x, —y,y), and 
consequently S satisfies E{—k,—l) = E{kJ). As a result: 

A+ = -A_. (B13) 


Using (B13) and shifting the axes in the resulting integrals 
{k —i' k + n/2 and 1^1 + ml2), reduces (Bll) to the 
following dispersion relation: 




dkdl K^{K^ 


N‘^)S^{k,l)x 


{mk — nl) [nm(A:^ — Z^) -|- (m^ — n'^)k+l+\ 
i/d (kK^ -{k + n)K‘^) + {a + 2r)K^K^ 
= 7r(cr -I- r) — iirnP, 


(B14) 


where A/ = {k + n)^ + {I + m)^. The corresponding disper¬ 
sion relation on a periodic box, can be readily calculated 
by simply substituting the integrals in (B14) by finite sums 
of integer wavenumbers. For a mirror symmetric forcing 
obeying: 

E{-kJ) = EikJ), (B15) 

the eigenvalues tr satisfy the symmetries (17). In order to 
show this, we consider (B14) for a(^-n,m) and change the 
sign of k in the integral to obtain: 




dkdl K^{K^ 


N‘^)S^{-kJ)x 


{mk — nl) [nm(A:^ — Z^) -I- (m^ — n^)fc+Z+] 

-i$ {kK"^ - {k + n) A2) -t (cr(_„_m) -h 2r) A/ 

= 'K{(T(-n,m) "I- r)N'^ + iTm/d. (B16) 


Taking the conjugate of (B16) and using the mirror sym¬ 
metry (B15) yields (B14) and therefore cr(_ri,m) = 
Similarly, it can be readily shown by considering (B14) for 
^(n,-m) and changing the sign of I in the integral, that 


in/3 — {a + r)N'^ = 

7V2 r°° r°° 

= / / [nm/k"^ — l'^) + {m'^ — n'^)kl] Snmdkdl 

= A+-A_, (Bll) 

where 

-| poo poo 

A± = — / / dkdlKl/Kl-N^)S^{k±,l±)x 

271" 7-00 7-00 

[nm(A:^ — P) + {mP — n^)ZcZ] {mk± — nl±) 

^ il3{k_Kl-k+Kp) + {a + 2r)KlKp ' 

(B12) 
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